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A method is presented for automatically calculating true first order rate constants for gas phase and wall 
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1. Introduction 

Steady-state tubular flow reactors are widely used in 
chemical kinetics studies as a means of converting the 
reaction time into a distance measurement. One way of 
using this reactor is to inject one reactant into a flowing 
carrier gas which contains a second reactant whose concen- 
tration is much greater than that of the first. The concentra- 
tion of the first reactant is then measured as a function of 
distance down the tube. If distance is assumed to be equal 
to the product of the reaction time and the average linear 
flow velocity of the carrier gas, then the concentration of the 
first species will be simply 



C - C e 



(1) 



where k is the first-order rate constant, (u) is the average 
carrier flow velocity, and Co is the value of C when the 
distance z is zero. The position used for the z origin is 
arbitrary since (1) shows that k can be determined from the 
relative concentration of C. However, it must be fixed at 
some distance downstream from the injection point so that 
the measurements start only after the reactants are well 
mixed. 

Unfortunately (1) is only approximate. Laminar flow exists 
in most experiments so that the carrier velocity will have a 
parabolic profile. Reactants near the tube center will travel 
faster than those near the wall. This will create a radial 
gradient in the concentration of C. The extent of this gradient 
is a complex function of the flow velocity, reaction rate, and 
molecular diffusion effects. An additional complication will 
arise if C can also be destroyed by a reaction on the tube 
wall. This is a common occurence if C is an atomic species. 
In spite of these complications, C still exhibits an exponential 
decay along the tube, as long as the observations are made 
sufficiently far downstream from the mixing region. Instead 
of (1), we now have 



C(r) = C (r)e- fc *' /<u> 



(downstream from the mixing region) 



(2) 



Here, the observed decay parameter A;* is a function of (u), 
the diffusion coefficient D c of species C in the carrier gas, 
the true first-order rate constant A;, and a rate constant k w for 



a first-order wall reaction; r is the distance from the tube 
center. This equation will hold at all points sufficiently far 
downstream from the mixing region. How far downstream is 
sufficient will be discussed later. The concentration in this 
region exhibits a radial distribution which remains constant 
along the tube. To determine k* correctly, it is of course 
necessary to measure C over the same range of r values at 
each particular position along the tube. 

Once A:* has been determined it must be related to k, the 
actual first-order rate constant. This problem has been 
solved by Walker [l]. 1 Unfortunately, the value of A:* is 
given by the first positive root of a polynomial having a 
fairly large number of non-negligible terms. Thus, it is not 
possible to give a closed expression for k or A* in terms of 
the other parameters. Walker has tabulated a few arrays of 
A:* values corresponding to a number of experimental condi- 
tions. To use his results, however, to correct experimental 
k* values is tedious since interpolation is required. Further- 
more, the accuracy of the interpolation has not been verified. 
The purpose of the present work is to make his method 
easier to use. To do this a simple computer program has 
been written in the form of a FORTRAN subroutine called 
ROOT which will calculate k or k* when given the values of 
the other parameters. A number of plots are presented of A:* 
versus k from which other values can be obtained by a linear 
interpolation whose accuracy has been verified. The program 
can also be used to determine higher order decay terms. 
From these it is possible to estimate the extent of the mixing 
region. 

Actually, Walker's solution need be used only if a wall 
reaction is present. When k w can be neglected, there exists 
[2-4] a very simple, and excellent approximate formula for 
k*. It is 



X* = V2(V(1 + 4X) - 1) 



(3) 



where X* = Gk*/(u)\ X - Gk/(u)\ and G = D c + a 2 (u) 2 / 
48D C ; a is the radius of the reactor tube. Later, a brief 
discussion regarding the origin of this formula will be given. 
A simplified derivation of Walker's solution is presented 
next. This will provide the basis for the discussion of the 
computer program. Following that an example demonstrating 
its use will be presented. 



1 Figures in brackets indicate the literature references at the end of this paper. 
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2. Solution of Diffusion Equation With 
First-Order Kinetics 

The partial differential equation which describes the 
tubular reactor with first-order kinetics under laminar flow 
conditions in the steady-state is 



2<u)(l - r 2 /a 2 ) 



dC 
dz 



= D r 



d 2 C IdC d 2 C 

di 2 rdr dz 2 



kC. 



(4) 



The first-order rate constant k w for the wall reaction appears 
in the boundary condition, 



/>, 



^u^r=a' 



(5) 



From elementary kinetic theory, k w = 1 /4yc, where y is the 
fraction of molecules of C which are destroyed on striking 
the surface, and c is the mean molecular speed of species C . 
(A more accurate formula for k w is given in ref. 5.) 

With the following dimensionless parameters, R — r/a, Z 
= z/a, D = D c /2a(u>, K = ak/(u), K w = k w /(u), (4) and (5) 
become 



(1 - R 2 



dC 
dZ 



D 



D 



d 2 C IdC d 2 C 
dR 2 + RdR + dZ 2 



1 /2K w C R=l . 



1 liKC (6) 



(7) 



In terms of these parameters, Walker's dimensionless param- 
eters are r* = r?, z* = Z, u = V2D, B 2 = K/2D, 8 
= 2D/K W 

The solution of (6) for the region downstream from the 
mixing point is 



C = 2 AigiWe-* 2 . 



(8) 



Values of the decay parameters Kf are fixed by the boundary 
condition (7). The appropriate mixture of the radial functions 
gi(R) could be determined if the radial concentration profile 
at the mixing point were known. However, if the measure- 
ment of C is begun at az value such that exp (—Kf) » exp 
( — Kf), for i ^ 2, then only the first term in (8) will be 
significant. The trne first-order rate parameter £ can thus be 
determined from the observed decay parameter Kf provided 
the functional relationship between the two can be deter- 
mined. 

To establish this relationship, begin by equating (2) to the 
first term in (8). This gives C(r) = C (r) exp (—k*z/(u)) = 
A lgl (R) exp (~KfZ) = A g(R) exp (-K*Z), which defines 
the dimensionless decay parameter K* = k*a/(u). Substitut- 
ing this into (6) and (7) yields a differential equation and 
boundary condition which must be satisfied by g(R). These 
are, 



#g 

dR 2 



+ 



RdR 



K* 2 + (1 - R 2 )K*/D - K/2D 

KwgR=\ 



2D 



g = (9) 



(10) 



This function can be expressed as a power series in even 
powers of R. 



g(R) = 2 B„R 2n . 



(11) 



The coefficients B n are given by 
B, = - l UaB 



B 



1 



" {2nf 



(K*/D)B 



„_ 2 - «#„_, 



(12) 



where 



a = K* 2 + K*ID - K/2D. 



When R = 0, we have g = B , so that B is proportional to 
the concentration at the tube axis. Since only the relative 
concentration is of interest here, B will be arbitrarily given 
the value unity. To derive (12) simply insert (11) into (9) 
and write out a few terms having different values of n 
starting with n = 0, and set the coefficients of like powers 
of R equal to zero. (If the more general power series 
containing both even and odd powers of R is used instead of 
(11), this procedure will show that only even powers of R 
have non-zero coefficients.) 

Equation (11) must also satisfy the boundary condition 
(10). Substituting it into (10) gives 



F(K* ,A:,A: W .,D) = 2 B n (2n + Kj2D) = F(x) = 0. 



(13) 



F is a function of £*, K, K w , and D. If the latter three 
quantities are specified, then the K* value of interest will be 
the first positive root of F considered to be a function of /C*. 
Alternatively, if £*, K w , and D are given, then K will be a 
positive root of F, but not necessarily the smallest one. In 
experiments where C is destroyed on the wall, K w can be 
determined by measuring K* in the absence of the second 
reactant. Then K = and (13) can be solved directly for£ M , 
in terms of K* and D to give 



K w — 



-2D 2 B n (2n) 

n=0 

00 



(14) 



The subroutine ROOT uses the Newton-Raphson method 
[6] to determine the roots of (13). If Xj is an approximate 



value of#, then a better value Xj +] is given by the recurrence 
relation 



«i+i = x j ~ F(xj)/F r (xj) 



(15) 



where F'(xj) is the derivative of F with respect to x at the 
valuer = Xj. The sequence x , Xj, . . . will converge to the 
correct root if the initial value is reasonably close to the 
correct value. 

The subroutine ROOT, shown in the appendix, requires 
the following calling statement, 
CALLROOT(ZS,Z,ZW,D,IOPT,F,B,NB,IFLAG) 



whe 



ZS = K*, the observed first order decay parameter, 
Z = A, the true first-order gas phase rate constant, 

ZW = Ky^ the first-order wall decay constant. 
D = D, the diffusion coefficient of species C in the 

carrier gas, 
IOPT - 1,2,3,4 



If IOFr = 1, ROOT evaluates the function F(K*,K,K W ,D) 
If IOPT = 2, then K* is calculated starting from an 

approximate value. 
If IOPT = 3, K is calculated from an approximate 

value. 
If IOPT = 4, K w is calculated. 
F = F{K*X,K W ,D\ defined by (13), 

B = B n the coefficients defined by (12). It is an array 
which must be dimensioned to 30 in the calling program. 



NB = the number of terms in (13) or (14) used by ROOT 
for a particular calculation. 

I FLAG = 1 if the number of iterations in the Newton- 
Raphson procedure exceeds the number IMAX. IFLAG = 
otherwise. 

To implement the Newton-Raphson method, ROOT cal- 
culates the derivative of F from the formula 



F = 2 #»(2n + KJ2D). 



(16) 



For IOPT = 2, we want the derivative with respect to A^*; 
here, the B n are given by 



B\ = -V 4 a' 

B 'n = ^-r 2 (B n -JD + K*B' n . s /D - a'fi n _, (17) 

a' = 2K* + l/D. 

For IOPT = 3, the derivative of F with respect to K is 
needed; in this case, the B' n are 



B\ = -Via' 
] 



B n — 



(IPffn-t/D-a'Bn^-aBn-r) 



(18) 



(2rc) 2 

a' = - V2D. 

ROOT continues to add terms to F until its absolute value 
changes by less than the number PREC1. Up to 30 terms 
are allowed by the dimension of the B array, but fewer than 
10 usually suffice for calculating the low order roots. (If the 
dimensions of B are changed in ROOT, it is also necessary 
to change the dimensions of the array BN which contains the 
derivatives of the B n values.) ROOT continues iterating until 
the absolute value of A* (or A) changes by less then the 
number PRFC2. A value of 1 X 10" 4 for both PREC1 and 
PREC2 has been used. The maximum number of terms used 
for calculating F is specified by the number NMAX, which 
must not exceed the dimensions oi B . The values of PREC1, 
PREC2, NMAX, and IMAX are assigned in the DATA 
statement. 

For ROOT to be automatic, i.e., for it to calculate A or 
A^*, given the other parameters, it is necessary to specify 
suitable starting values for the Newton-Raphson calculation 
in terms of these parameters. When A* is desired, (IOPT = 
2), a satisfactory initial value is 0.9 A. For a determination 
of K (IOPT = 3), the initial value 1.2 K* is satisfactory. 
These starting values have been checked and found to give 
convergence to the correct root over the following ranges of 
parameter values; 



^ K* ^ 0.68 
0^/C^ 0.58 
^ K w ^ 0.2 
0.01 ^D^ 1.0 



(19) 



They may also be satisfactory outside these ranges, but 
verification of this would require further investigation. 

To show the functional dependence of A^* on the other 
parameters, a series of plots of A* versus K is given in 
figures la through 6b. The plots are arranged in pairs. Each 
pair contains plots having the same value of K w . The first 
figure in the pair contains plots for different D values 
running from D = 0.01 to 0.07, and the second contains 
these where D = 0.10 to 1.00. The six values of K w were 
0.0, 0.01, 0.02, 0.06, 0.10, and 0.20. Values of K* 
corresponding to K w values lying between those given, can 
be obtained by a linear interpolation whose accuracy is 
better than 1%. This probably exceeds the accuracy with 
which A* can be determined from the figures. The ranges 
covered by the parameters in these plots are the same as 
those given by (19) and should be sufficiently large to 
encompass the conditions encountered in the majority of low 
pressure, high flow velocity experiments. 

3. An Example Illustrating the Use of ROOT 

To demonstrate the use of ROOT, consider the following 
experiment. There is the reaction A + B — > products, with a 
second order rate constant k 2 and with [A] « [B] so that a 



first order rate constant k = k 2 [B] can be defined. The 
reactants are contained in a carrier gas flowing through a 2 
cm I.D. tube with an average linear velocity of 200 cm/s. 
The measurement of the concentration of A is begun 5 cm 
downstream from its injection point and is found at 20 cm to 
have decayed exponentially to 0. 1 of its value at 5 cm. 
When species B is removed from the carrier, A is observed 
at 20 cm to have decayed to 0.3 of its value at 5 cm. Since 
A decays in the absence of B, it is assumed to be undergoing 
a first order wall reaction. Thus, k* (total) = 30.7 s _1 , and 
k* (wall) = 16.05 s _1 . If the diffusion coefficient of A in the 
carrier gas at the prevailing pressure is 20 cm 2 s _1 , then 
the values of the dimensionless parameters used by ROOT 
are, K* (total) = 0.1535, K* (wall) = 0.0803, and D = 
0.05. 

The first task is to use ROOT to determine K w from the 
value of K* observed in the absence of B . The following 
calling program will accomplish this. 

PROGRAM CALL 

DIMENSION B(30) 

ZS = 0.0803 

Z= 0.0 

D = 0.05 

I0PT= 4 

CALL ROOT(ZS,Z,ZW,D,IOPT,F,B,NB,IFLAG) 

PRINT ZW 

END 

In this case ROOT will return a value of 0.05 for ZW 

To determine K from K* (total), use K w = 0.05, and use 
1.2 K* = 0.1842 as an approximate value of K. Then use 
ROOT with IOPT = 3 to determine the correct value of K. 



The following program accomplishes this. 

PROGRAM CALL 

DIMENSION B(30) 

ZS = 0.1535 

Z = 1.2*ZS (approximate value) 

ZW = 0.05 

D = 0.05 

IOPT = 3 

CALL ROOT(ZS,Z,ZW,D,IOPT,F,B,NB,IFLAG) 

PRINT Z 

END 

For this example, ROOT gives 0.0842 for the accurate value 
ofK. 

From these results, the values, k = 16.84 s , and k w = 
10.00 cm s _1 are obtained for the laboratory parameters. 

It has been assumed in this experiment that the reactants 
are well mixed 5 cm downstream from the injection point of 
A. This is equivalent to assuming that only the leading term 
in (8) is being observed. ROOT can be used to determine 
the higher order decay parameters in (8). To do this for the 
above example, the values K = 0.0842, K w = 0.05, and D 
— 0.05 are used with IOPT = 1 to calculate F as a function 
of K*. From this the zero crossing points which correspond 
to higher roots of (13) can be located and used with IOPT = 
2 to get accurate values. For the present example, the next 
decay parameter is K$ = 1.3144, which yields k* = 263 
s _1 . Therefore, at 5 cm downstream from the mixing point, 
the second term in (8) would have decayed to 0.0014 of its 
initial value. This observation does not guarantee that only 
the leading term in (8) is being observed. If the coefficient 
A 2 were much larger than A x then the second term in (8) 
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Figure la. Plots ofK* versus Kfor D = 0.01, 0.02, 0.03, 0.04, 0.05, Figure lb. Plots of K* versus Kfor D = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 

0.06, 0.07. 0.7, 0.8, 0.9, 1.0. 

K u ,= 0.0. K w =0.0. 
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Figure 2a. Plots of K* versus Kfor D = 0.01, 0.02, 0.03, 0.04, 0.05, Figure 2b. Plots of K* versus K for I) = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 

0.06, 0.07. 0.7, 0.8, 0.9. 1.0. 

K w =O.Ol. K u = ().0\. 
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Figure 3a. Plots of K* versus Kfor D = 0.01, 0.02, 0.03, 0.04, 0.05, Figure 3b. Plots of K* versus Kfor D =0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 

0.06, 0.07. 0.7, 0.8, 0.9, 1.0. 

K W =()A)2. K w =0.02. 



could still be significant even though its exponential factor 
were small. As mentioned earlier, the coefficients A * could 
be determined from the radial distriubtion of C at the mixing 
point. This distriubtion is, however, rarely known. Pirkle 
and Sigillito [7] have solved (4) for several cases in which 
the initial distribution of C was uniform. In none of them did 
the values of A*, i = 2, exceed those oiA x . Unfortunately, it 



is difficult to say what the situation would be for other initial 
distributions. In practice, a logarithmic plot of the concentra- 
tion versus distance is usually made and only the linear 
portion (within experimental error) is used to determine A;*. 
Such a procedure, coupled with an examination of the next 
higher decay parameters k$ is most likely adequate to 
ensure observation of only the leading term in (8). 
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Figure 4a. Plots of K* versus Kfor D = 0.01, 0.02, 0.03, 0.04, 0.05, 
0.06, 0.07. 

A' u .= 0.06 



Figure 4b. Plots of K* versus Kfor D = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 
0.7, 0.8, 0.9, 1.0. 

K w = 0.06. 



























































D 


























,/ 


& 


0.07 










K w = 


0.1 












■it- 


^ 


/ 






















* 


4 


V\ 


S s 


/ 




















,<£ 


# 


y A 


y 


y 


s 


0.01 




K* 












^ 


% 


-/ 
























4 


P 


^ 


s 




^ 
















S< 


4 


^ 


y 




y 


















y/y 


P 


y 


y 






















d 


5 


y y 




/ 






















v 


y 


y 


























/ 




























0- 




























(a) 



























































D 
0.1 




0.60- 






















A 


y 












K w 


« 0.1 












r / 


Y 


r 


























/, 


2 


-**' 




















y 


l 


^ 









1.0 




K* 










.tt 


4 


^ 


^ 


.£> 




















A 


4 


i 





<£^ 




















S& 


%yy 


yy?' 
























A 


w 


'^y 


























w 






















































































0- 




























(b) 



0.20 



0.40 



0.60 



0.20 



0.40 



0.60 



Figure 5a. Plots ofK* versus Kfor D = 0.01, 0.02, 0.03, 0.04, 0.05, Figure 5b. Plots of K* versus Kfor D = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 

0.06, 0.07. 0.7, 0.8, 0.9, 1.0. 

a:„,= o.io. k mj =o.io. 
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Figure 6a. Plots of K* versus K for I) = 0.01, 0.02, 0.03, 0.04, 0.05, Figurk 6b. Plots of k* versus K for I) = O.l, 0.2. 0.3, 0.4, 0.5, 0.6, 

0.06, 0.07. 0.7, 0.8, 0.9, 1.0. 

K u . = 0.20. A.,,. = 0.20. 



4. Determination of k in the Absence of a 
Wall Reaction 

In the absence of any chemical reactions, Taylor [2, 3] 
and Aris [4] have shown that the coupling of axial and radial 
molecular diffusion with the parabolic velocity profile can 
be described by using the pseudo-diffusion coefficient G, 
given following (3), and considering only axial diffusion. 
The physical basis for the functional dependence of G on the 
molecular diffusion coefficient D c can be seen by considering 
a species initially distributed in a radial plane moving with 
the average velocity of the carrier gas. As this distribution 
moves down the tube it will at first take the shape of a 
paraboloid. Molecules near the tube wall lag behind those at 
the center. However, molecular diffusion in the radial 
direction allows these molecules to move away from the wall 
into the faster flowing carrier. There will therefore be a net 
movement of the distribution down the tube. For small 
values of D c , the distribution becomes dispersed about the 
radial plane just as though it were undergoing axial diffusion 
only, with a diffusion coefficient given by the second term in 
G. As D c becomes larger, the axial dispersion will decrease 
until a point is reached where true molecular axial diffusion 
surpasses the apparent axial diffusion resulting from the 
interaction of radial molecular diffusion with the velocity 
profile. Then the first term in G becomes dominant. This is 
simply Z) c , the molecular diffusion coefficient. 

When the diffusing species is being destroyed by a first 
order gas phase reaction, it is a good approximation [8] to 
consider only axial diffusion and use G for the diffusion 
coefficient. This yields the simple differential equation 



dz dz 2 



kC. 



(20) 



For the boundary conditions, C(z = 0) = Co, and C(z — > °°) 
— » 0, the solution of this equation may be written as 



c = C„c 



-k*zl<u> 



(21) 



To establish the relationship between A* and A\ evaluate the 
derivatives of C from (21) and substitute them along with 
(21) into (20). This yields the quadratic equation 

A* 2 + ({uf/G)k* - ((uf/G)k = 0. (22) 

Solving for A;* gives 

h* = l /2(-<u>VG ± V {({uf/Gf + 4((uf/G)k}). (23) 

To get (23) into the form of (3), multiple it by G/(u) 2 . This 
yields 

Gk*/(u) 2 = V 2 (-l ± V{1 + 4Gk/(u) 2 }) = X*. (24) 

Since X* must be positive to satisfy the boundary condition 
C(z — > °°) — > 0, the positive value of the square root must be 
used. The resulting expression is (3). Values of A* given by 
(3) were compared to the exact ones determined by Walker's 
method for the range of parameter values shown in figures 
la and lb. The agreement averaged better than 0.2 percent 
over the whole range. 



5. Appendix. Listing of Subroutine ROOT 

The program is shown here in single preeision. To calcu- 
late decay parameters higher than K$ , rewrite it in double 
preeision and increase the dimensions of the arrays B and 
BN to 40. 



10 

11 
12 
13 
14 
15 
IS 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
38 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 



SU6EQUTINE RO0TCZS.Z.ZU.D. I0PT.F.6.NB. IFLPG) 
DIMENSION BC30).DBC30) 

DATA PREC1.PREC2.NMAX. IMAX/.0001. .0001.30. 10/ 
IFLAG=0 
R2D-i./(2.*B) 
ZUD=ZU*R2D 
RD-l./D 
I TAB =9 
70 ZSD=ZS*RD 

A=ZS*ZS+ZSD-Z*R2D 

IFCIGPT.EQ.2) DA»2.*ZS+RD 

IF': l OPT. EQ. 3) DA— .5#RD 

6< !)=--. 25*A 

DB(D=-,25*DA 

B<2) ". 9S25*(ZSD-A*B CD) 

IF(iOPT.EQ.2) DB(2)=.0S25*(RD-DA*BC1)-A*DB(1)) 

IFCI0PT.EQ.3) DB(2)=.0625*(-DA*B(1)-A*DB<1)) 

F-ZUD 

DF = ,0 

GN-.0 

GD = i. 

DO 15 N = 1.2 

D 1 =2*N 

Q2=Q1-ZUD 

IFU0PT.EQ.4) GO TO IS 

F»F+B(N)*G2 

IF(CI0PT.EQ.2) .OR. (I0PT.EQ.3)) DF =DF+DB(N)*G2 

GO TO 15 
16 GD=GD+B(N) 

GN=GH+B(N)*G1 
15 CONTINUE 

DO 10 N*3,NMAX 

NB=N 

Q 1 -2*N 

Q2»Q1*9! 

RG2=1./Q2 

B(N)»RQ2*(ZSD*BCN-2)-A*BCN-l)) 

GO TO (11. 12. 13.. 20). I0PT 
12 DB CN) =RQ2*(RD*BCN-2)+ZSD*DBCN-2)-DA*B (N- 1 ) -A*DB (N- 1 ) ) 

GO TO 11 



41 


13 


jJbi.N;=Ra2*(ZSD*DB(N-2)-Dft*3(N-l)-A*DB(N-l)) 


42 


11 


Q3=Qi+ZUD 


43 




FOLD-F 


44 




F=F+B(N)*Q3 


45 




IFCCIGPT.EQ.2) .OR. ( I0PT.EQ.3) ) DF=DF+DB(N)*Q3 


4S 




IF(ftBSCF-FOLD).LT.PRECl) GO TO 21 


4? 




GO TO 19 


48 


20 


GD=GD-i-BCN) 


49 




GN0LD=GN 


53 




GN=GN+8(N)*Q1 


5i 




IF(ABS(GN-GNOLD).LT.PRECn GO TO 21 


52 


10 


CONTINUE 


53 


21 


GO TG (39.30.32,56), IOPT 


54 


30 


ZSGLD=ZS 


55 




ZS--ZS-F/DF 


56 




IF(hBS(ZS-Z50LD).LT.PREC2) GO TO 99 


57 




GO TO 31 


53 


32 


ZOLD-Z 


59 




Z=Z-F/DF 


60 




IF(ABS<Z-Z0LD).LT.PREC2) GO TO 99 


61 


31 


ITAB=ITAB+1 


62 




IFCITAB.GT.IMAX) GO TO 41 


63 




GO TO 70 


64 


50 


ZU»-2.*D*GN/GD 


65 




GC TO 99 


66 


41 


IFLAG=j 


67 


99 


RETURN 


-P 




END 
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